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INTRODUCTION 


The  R-matrix  propagation  technique  has  recently  been  introduced  by  several  authors 
(References  1  through  6)  as  a  means  of  removing  the  numerical  instability  encountered  in 

calculating  diffi^ction  from  deep  gratings  (depth  h  >  wavelength  A,).  The  method  was  first 
described  in  the  field  of  chemical  physics  (Reference  7).  The  idea  is  to  find  a  quantity 
called  R  that  describes  the  relationship  between  some  field  quantity  and  its  derivative. 

When  the  field  has  exponential  behavior,  exp(az),  the  quantity  R  will  be  proportional  to  az 
and  will  not  suffer  fi’om  exponential  instability.  Applying  this  idea  to  such  a  relationship  in 
the  theory  of  grating  diffraction,  the  method  seeks  a  matrix  that  relates  the  electric  field  on 
both  sides  of  a  layer  to  the  magnetic  field  on  both  sides.  This  method  on  the  R-matrix 
method  is  unlike  the  T-matrix  propagation  method  that  seeks  a  matrix  relating  the  electric 
and  magnetic  field  on  one  side  of  a  layer  to  the  electric  and  magnetic  field  on  the  other  side 
of  a  layer. 

In  References  1  and  2,  diffraction  calculations  were  performed  for  nonabsorbing 
dielectric  and  absorbing  metallic  sinusoidal  gratings  with  period  d  =  1.7A,  and  heights 

h  =  0.17,  1.7,  17,  and  170X  for  both  transverse-electric  (TE)  and  transverse-magnetic 
(TM)  polarization.  The  results  were  compared  with  other  methods  where  apphcable.  The 
performance  of  the  R-matrix  method  was  good  in  all  cases  except  for  TM  polarization  and 

metallic  gratings  with  h>  17X.  These  cases  suffered  from  poor  convergence  where  the 
diffraction  efficiency  fluctuated  with  the  number  of  orders  kept  in  the  expansion. 

In  Reference  1,  the  R-matrix  method  was  used  in  conjunction  with  a  multilayer-modal 
expansion  to  calculate  the  diffraction  efficiency.  The  grating  is  subdivided  into  many 
sublayers  such  that  within  each  sublayer,  the  permittivity  can  be  approximated  as  constant 
across  the  sublayer.  In  the  dimensions  parallel  to  the  boundaries  of  a  sublayer,  the 
permittivity  and  the  fields  are  expanded  in  a  Fourier  series.  Within  each  sublayer,  the 
Fourier  components  of  the  permittivity  are  known,  and  the  Fourier  components  of  the  field 
are  determined  from  a  set  of  linear  equations. 

In  this  work,  we  do  not  expand  the  permittivity  and  fields  in  each  sublayer  into  Fourier 
components.  Instead,  we  will  solve  for  the  fields  in  real  space.  The  motivation  for  this 
approach  is  to  try  to  circumvent  the  slow  convergence  of  a  Fourier  expansion  when  there 
are  one  or  more  sharp  discontinuities  in  the  quantity  being  expanded.  TTiese  discontinuities 
include  the  field  across  a  boundary  or  the  sharp  jump  in  the  dielectric  permittivity.  It  is  not 
clear  to  the  authors  why  metallic  permittivities  and  TM  polarization  appear  to  present  the 
greatest  problem  in  this  regard.  Nevertheless,  this  problem  has  motivated  taking  the 
present  real  space  approach  where  series  representation  of  discontinuities  in  permittivity  is 
not  an  issue. 
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The  method  described  here  is  similar  to  that  described  in  References  3  and  4  where  the 
authors  mainly  looked  at  transmission  and  band  stracture  of  a  bulk  photonic  crystal  and 
dispersion  of  surface  modes  propagating  along  the  surface  of  a  truncated  photonic  crystal. 
These  calculations  were  for  nonabsorbing  dielectric  materials.  In  this  paper,  we  will  focus 
on  the  gratings  considered  in  References  1  and  2.  The  next  section  briefly  describes  the 
method  used  in  this  work.  Additional  details  can  be  found  in  References  1,  3,  and  4. 
Numerical  results  along  with  convergence  data  are  presented. 


METHOD  OF  CALCULATION 


In  Figure  1,  a  sinusoidal  profile  with  nomenclature  is  shown.  We  assume  the  grating 
profile  is  infinitely  periodic  in  the  x  direction,  uniform  in  the  y  direction,  and  with  peak- 
to-valley  height  /z  =  Z2  -  Zj  in  the  z  direction.  The  grating  is  illuminated  by  a  plane  wave 
of  wavelengdi  A  with  wave  vector  in  the  (x,z)  plane.  There  are  three  general  regions  of 
interest:  the  homogeneous  superstate,  the  homogeneous  substrate,  and  the  finite  thickness 
region  in  between  containing  the  grating  profile.  The  superstate  has  permittivity  and 
contains  the  incident  and  reflected  fields.  The  substrate  has  permittivity  and  contains 
the  transmitted  field.  The  profile  region  in  between  is  described  by  a  spatially  variable 
permittivity  £(r)  where  r  =  ix,z)  and  Zj  <  z  <  Z2. 

In  Figure  1,  a  sinusoidal  profile  is  subdivided  into  layers  each  having  uniform 
thickness  Az  =  hi n^.  Even  though  the  grating  profile  region  is  generally  divided  into 
sublayers,  the  number  of  sublayers  can  range  from  1  to  much  greater  than  1.  For  a 
given  height  h,  the  number  obviously  relates  to  thickness  Az,  and  the  main  criterion  for 
Az  is  the  requirement  that  the  permittivity  within  each  sublayer  can  be  treated  as 
independent  of  z.  For  example,  a  rectangular  profile  need  not  be  subdivided 
regardless  of  height  h.  However,  for  a  sinusoidal  profile,  it  is  normally  necessary  to 
subdivide  into  /J2»1  layers  to  achieve  a  reasonably  accurate  profile  shape  and  sublayers 
which  are  approximately  independent  of  z.  Even  though  we  assume  each  sublayer  has 
constant  thickness  Az,  this  assumption  is  not  necessary. 
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Superstrate 


Substrate 


FIGURE  1.  Schematic  of  Sinusoidal  Profile  Versus  One  Period  Showing  x- 
Coordinate  Discretization  and  Sublayer  Divisions.  In  this  example,  the  dots  indicate 
discretization  of  the  j:-dimension  over  one  period  d  into  =  29  points.  The  number 
of  divisions  of  grating  height  his  /Zz  =  10,  as  shown  by  the  horizontal  solid  lines. 
The  area  under  the  sine  curve  contains  substrate  materi^.  Each  dot  lying  under  the 
sine  curve  represents  an  jr-point  where  the  permittivity  is  that  of  the  substrate. 
Likewise,  each  dot  above  the  sine  curve  is  an  jf-point  with  superstrate  permittivity. 
Each  adjacent  pair  of  horizontal  soMd  lines  represents  a  sublayer,  and  the  dotted  line 
in  between  also  denotes  the  center-z  coordinate  for  that  sublayer. 


We  wish  to  solve  the  two  Maxwell's  equations  VxE(r)  =  z(m/c)H(r)  and 
VxH(r)  = -i(ta/c)s(r)E(r)  within  a  given  sublayer.  To  do  this,  we  write  the  x- 
dependence  of  these  two  equations  as  a  centered  finite-difference  approximation.  The  x- 
coordinate  is  discretized  over  period  d  into  points  each  uniformly  separated  by 
/!x)c  =  dln^.  For  a  layer  bounded  from  z  — >  z  -l-  Az ,  we  eliminate  the  z-component  of  the 
electric  and  magnetic  fields  and  obtain  four  equations  relating  the  x  and  y  field 
components: 


dE^{x,z)  _  ico 


dz 


Hy  (_X,  z) 


ic  {Hy(x  +  Ax,z)-Hy(x,z)  Hy{x-Ax,z)-Hy{x,z)' 

6)(Ax)^  1  £(jc  +  Ax/2,z^)  e(x-Ax/2,z^) 


(la) 
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dEy{x,z)  i(o 

— - = - 

oz  c 


dHx(x,z)  _  id) 


= - e(x,z^)Ey(x,z) 

dz  c  ^ 

^Ey  (x,z)-  Ey  (x  +  Ax,  z)-  Ey  (x  -  Ax,  z)  j 


(0(Ax) 


(lb) 


(Ic) 


dHy(x,z) 

dz 


1(0 


£{X,Zr)EJx,z) 


(Id) 


In  Equations  la  and  Ic,  the  finite-difference  centered  derivative  approximations  on  the 
right-hand  sides  are  accurate  to  0(Ax^).  As  shown  in  the  Appendix,  we  have  also 
considered  approximations,  which  are  accurate  to  0(Ax‘^)  (shown  in  the  Appendix).  We 
have  evaluated  the  permittivity  e  at  z  =  z^,,  which  represents  the  center  z-coordinate  of  the 

sublayer.  Since  the  sublayer  is  bounded  from  z  to  z+Az,  it  follows  that  z^  =z  +  Az/2. 
Because  of  this,  these  coupled  differential  equations  have  coefficients  which  are 
independent  of  z  for  a  given  sublayer.  When  all  discrete  x-coordinates  (denoted  by  X) 
are  included,  the  coupled  set  of  differential  equations  in  Equation  1  may  be  written  in 
matrix  form  as 


f£jx,zy 

r^.(X,z)^ 

Ey(X,Z) 

Ey(X,z) 

H,iX,z) 

=  M(X,z,) 

HAX,z) 

Hy(X,Z)y 

Hy(X,Z)y 

(2) 


Since  M  is  independent  of  z  within  the  sublayer,  the  modal  solution  is  straightforward  by 
diagonalization  of  M  and,  in  the  original  basis  set,  this  solution  has  the  form 


(E,(X,z)^ 

Ol2 

Ey(X,z) 

H,(X,z) 

= 

•^21 

exp(AiZ)Ci  -t- 

■^22 

exp(A2z)C2  -!-••• -i- 

Ein 

exp(Ayvz)C;v  (3a) 

^Hy(X,z)j 

[EmJ 

1*5^2  J 
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where  N  =  4n^ .  The  column  vectors  on  the  right-hand  side  are  eigenvectors  of  the  matrix 
M(X,z^),  and  the  Xj,  j  =  1,  are  the  associated  eigenvalues.  The  Cj  are  constants. 

Equation  3a  may  be  compactly  written  in  matrix  form  as 


Ey{X,z) 

H,(X,z) 


S(X,zJexp(Az)C 


(3b) 


The  S(X,Zc)  is  a  square  matrix  having  columns  which  are  the  eigenvectors  of  M(X,z^). 

The  exp(Az)  is  a  diagonal  matrix  of  exponential  terms  with  A  representing  the  set  of 
eigenvalues  associated  with  M(X,z^),  and  C  is  a  column  vector  of  constants  Cj. 

In  constructing  the  matrix  M,  field  terms  in  Equations  la  and  Ic  with  extreme  x- values 
(xl  Ax),  which  fall  outside  the  dimension  of  one  period  are  "wrapped  around"  by  using 
the  Floquet  relationship.  Because  the  only  part  of  Equation  1  that  depends  on  the  grating 
profile  is  the  permittivity  £(r),  the  task  of  describing  various  profile  shapes  is 
straightforward,  and  typically  only  minor  changes  to  computer  algorithms  are  required. 

From  Equation  3b,  we  may  form  a  relationship  which  is  the  essence  of  the  R-matrix 
algorithm 


E.a^^z)  ^ 
Ey(X,z) 
E^(X,z  +  Az) 


( 


=  r(Az) 


H,(X,z)  ^ 
Hy{X,z) 
H,{X,z  +  ^z) 


[Ey(X,z  +  Az)) 


[Hy(X,z  +  Az)) 


(4) 


This  is  accomplished  by  combining  Equation  3b  evaluated  at  z- values  of  z  (as  shown)  and 
z  +  Az  to  eliminate  C  and  rearranging  the  result  to  have  the  form  given  in  Equation  4.  This 
matrix  equation  defines  the  sublayer  r-matrix  and  is  basic  to  the  R-matrix  algorithm.  The  r- 
matrix  must  be  calculated  for  each  sublayer.  We  compare  the  r-matrix  definition  to  the 
corresponding  definition  for  a  sublayer  t-matrix  which  we  write  as 


^E,{X,z)'^ 

f 

Ey(X,z) 

H,(X,z) 

=  t{Az) 

^Hy(X,Z)y 

\ 

£/X,z  +  Az) 
H^(X,z  +  Az) 
HJX,z  +  Az) 


(5) 
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The  elements  of  the  t-matrix  contain  terms  which  behave  exponentially  as  exp(±aAz) 
where  a  is  complex,  and  Az  is  the  thickness  of  the  sublayer.  When  the  real  part  of  the 
product  aAz  is  large,  numerical  problems  with  exponential  overflow  and  underflow  will 

likely  occur.  Even  if  the  real  part  of  aAz  is  small  for  any  sublayer,  propagation  through 
many  sublayers  by  means  of  products  of  t-matrices  will  cause  the  exponential  behavior  to 
accumulate  with  the  number  of  sublayers,  which  will  eventually  lead  to  numerical 
instability.  On  the  other  hand,  the  elements  of  the  r-matrix  have  a  linear  dependence  on 

±aAz,  and  this  behavior,  which  is  a  vast  improvement  over  exponential,  will  be  much 
more  numerically  tractable  for  propagation  through  many  sublayers. 

Returning  to  Equation  4,  note  that  if  Az  is  the  full  thickness  h  of  the  profile  region  (for 
example,  a  rectangular  profile),  then  no  further  propagation  through  the  profile  region  is 

necessary  because  the  fields  have  been  related  by  r(Az)  across  /?  =  Z2  -Zi-  For  more 
general  profile  shapes,  more  sublayers  are  needed.  Then  it  is  necessary  to  relate  the  field 

solutions  obtained  in  each  sublayer  of  thickness  Az  that  subdivides  =  Z2  -  Zi-  To  this 
end,  we  also  assume  that  a  relation  analogous  to  Equation  4  exists  that  covers  the  entire 
grating  height  h  =  Zi-Z2-  We  express  this  relationship  as 


r^.(x,zi)^ 

r^.(x,zi)^ 

Ey(X,Z,) 

^Ri](Z2  Zj)  ^12(^2  ^1)'' 

H/X,z,) 

^.(X,Z2) 

vl^21  (^2  “  ^1 )  ®22  (^2  “  )> 

H,(X,Z2) 

^Ey{X,Z2)j 

^Hy(X,Z2)j 

where  the  global  R-matrix,  R(z2  -Zj),  is  written  in  a  sectored  form.  The  matrices  Rjj, 
Rj2,  R21,  and  R22  are  computed  by  means  of  a  recursive  algorithm,  and  further  details  are 
given  in  References  1  through  3.  To  initialize  the  recursion  relations,  we  see  from 
Equations  4  and  6  that  we  may  set  Z2  =  Zi  +  Azj ,  and  therefore  R( Azj )  =  r( Azj ) .  With  this 
initialization,  successive  applications  of  the  recursive  algorithm  yields  R(Azj -i- AZ2), 
R(Az,  +  Az2  +AZ3),  R(^i  +Az2  +  Az3  +-"  +  Az„P  where  Azj  +  Az2  +  AZ3  +  •  •  • -1- Az„^ 
=  Z2  -  Zi  and  the  Azj  can  be  different. 

Given  the  R-matrix  and  a  few  more  intermediate  steps,  we  can  arrive  at  the  final 
relationship  from  this  calculation  the  relative  amount  of  energy  diffracted  into  the  substrate 
and  superstate.  To  do  this,  we  convert  the  fields  to  ^-space  by  applying  a  Fourier 
transform  matrix  to  Equation  6.  This  step  is  advantageous  because  in  k-space  we  can  relate 
the  electric  and  magnetic  fields  in  the  homogeneous  substrate  and  superstate  (which  will  be 
important  later  in  Equation  9).  Another  advantage  is  that  grating  diffraction  is  conveniently 
expressed  in  k-space.  The  Fourier  transform  matrix  F(K,X)  is  discretized  consistent  with 
the  r-space  digitization  earlier  and  for  electric  fields  we  have 
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F(K,X) 


rF,(x,zi)^ 

rF,(K,Zi)^ 

E,(X,Zi) 

F/K,zi) 

E,(X,Z2) 

F,(K,Z2) 

^E^(X,Z2)^ 

,F,(K,Z2), 

(7) 


and  similarly  for  the  magnetic  field.  The  vector  K  denotes  the  set  of  wave  vectors  in 
^-space,  and  the  number  of  x-points  is  identical  to  the  number  of  diffracted  orders  in  the 
set  K.  We  apply  this  F  matrix  to  Equation  6,  which  yields 


fE,(K,Zrr 

Ey(K,zO 

F,(K,Z2) 


=  R(Z2-Zi) 


^H,(K,zO 

Hy{K,z,) 

H,{K,Z2) 

^Hy{K,Z2) 


\ 


(8) 


where  R  =  FRF“^ .  Since  the  Zj  and  Z2  are  the  substrate  and  superstate  boundaries, 
respectively,  we  can  use  the  boundary  conditions  to  relate  the  substrate  and  superstate 
fields.  From  Equation  8,  we  find 


^  F,'(K,Zi)  ^ 

f;(K,Zi) 

E/(K,Z2)  +  F;'’^(K'”^Z2) 
E/(K,Z2)  +  Fr(K“^Z2) 


V  3' 


=  R(Z2-Zi) 


H,'(K,Zi) 

77;(K,z,) 

H/(K,Z2)  +  7f/^(K'-,Z2) 
^H/(K,Z2)  +  /f;"^(K'"^Z2)^ 


(9) 


where  K"“  is  the  in-surface  incident  beam  wave  vector.  The  incident  beam  electric  and 
magnetic  fields  are  described  by  E'""  and  =(///",//;"").  The 

superstate  and  substrate  media  contain  the  reflected  and  tansmitted  fields,  denoted  by  r 
and  t,  respectively.  At  this  point,  a  relation  between  the  electric  and  magnetic  fields  can  be 
easily  obtained,  and  the  calculation  proceeds  exactly  as  described  previously  (Reference  1 , 
Equations  12  through  16). 


NUMERICAL  RESULTS 


The  numercial  results  of  this  paper  are  an  extension  of  Reference  1.  This  is  because  the 
numerical  results  of  this  work  are  also  based  on  the  same  grating  design  parameters 
considered  in  References  1  and  2  but  using  a  real-space  modal  expansion  method.  We  will 
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compare  some  results  of  the  present  work  with  previous  results  by  LP  and  the  authors' 
^-space  approach.  Figures  2  through  7  show  our  r-space  convergence  data  for  TM  and  TE 
polarization  (electric  field  parallel  and  perpendicular  to  the  plane  of  incidence)  for  three 
different  gratings.  The  figures  show  relative  intensity  versus  number  of  diffracted  orders. 
The  number  of  orders  is  equivalent  to  the  number  of  discretized  x-points  which  we 
choose  to  be  an  odd  number,  and  the  relative  intensity  refers  to  the  ratio  of  diffracted  power 
to  incident  power.  For  all  the  cases  considered,  the  number  of  sublayers  is  50,  and  the 
incident  angle  is  30  degrees. 

In  Figures  2  and  3,  convergence  data  is  shown  for  three  reflected  orders  of  an 
absorbing  metallic  grating  =  (-48.91,  4.2)  with  height  h  =  \.1X  and  period  d  =  1.1% 

{hld=\).  The  data  in  Figures  2  (TM  polarization)  and  3  (TE  polarization)  indicate  excellent 
convergence  versus  number  of  orders  for  this  grating.  For  the  same  grating  and  TM 

polarization,  Li^  (Reference  2)  reported  that  his  method  has  poor  convergence,  because  the 
results  keep  fluctuating  with  the  number  of  orders.  In  Table  1,  we  compare  the  average 
value  of  the  last  five  r-space  data  points  shown  in  Figures  2  and  3  with  previous  results 
given  in  References  1  and  2.  The  agreement  for  TE  polarization  is  veiy  good.  The 
agreement  for  TM  polarization  is  poor,  but  that  is  to  be  expected,  because  Ae  previous 
results  have  difficulty  with  convergence.  For  h/d  =  1,  the  present  results  compare  well 
with  the  integral  meAod  and  the  method  of  fictitious  sources  (Reference  8),  which  gives 
TM-polarized  reflection  diffraction  eeficiencies  for  the  -2,-1,  and  0  order  as  0.197,  0.086, 
and  0595,  respectively. 


FIGURE  2.  Relative  TM-Polarized  Diffraction  Intensity  of  the  Reflecting 

0,-1,  and  -2  Orders  Versus  Number  of  Orders.  The  grating  profile  is 
sinusoidal,  and  the  grating  material  is  an  absorbing  metal  with  permittivity 

(-48.91,  4.2).  The  grating  height  is  h  =  \.lX,  period  d  =  \.lX,  angle  of 
incidence  30  degrees,  and  the  number  of  profile  sublayers  =  50. 
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FIGURE  3.  Relative  TE-Polarized  Diffraction  Intensity  of  the 
Reflecting  0,-1,  and  -2  Orders  Versus  Number  of  Orders.  This  case 
is  identical  to  Figure  2  except  the  incident  beam  is  TE  polarized. 


TABLE  1.  Reflection  Diffraction  Efficiencies  Absorbing  Metallic  Sinusoidal  Grating. 


il 

3 

□ 

o 

11 

3 

Order 

5“point 

r-space 

average 

Our 

previous 

fc-space 

Li’s 

1 

5-point 

T'Space 

average 

Our 

previous 

^-space 

Li’s 

-2 

0.4122 

0.4152 

0.4135 

0.1914 

0.1885 

-1 

0.3348 

0.3338 

0.3353 

0.1365 

0.1353 

0 

0.2039 

0.2010 

0.2018 

0.3093 

0.3086 

■■nxjgi 

-2 

0.1560 

0.0274 

0.1264 

0.2834 

0.0160 

-1 

0.0705 

0.0664 

0.0603 

0.0462 

0.0078 

0 

0.5993 

0.2146 

0.6609 

0.1451 

0.0315 

0.1618 

In  Figures  4  and  5,  we  consider  the  same  absorbing  metallic  grating  except  that  the 
height  is  10  times  greater.  The  grating  height  h  is  now  17X  (h/d  =  10)  while  the 
permittivity  and  period  d  remain  unchanged.  Convergence  data  for  the  three  reflected 
orders  are  shown  in  Figures  4  (TM  polarization)  and  5  (TE  polarization).  Both  cases 
exhibit  a  "slow"  damped  oscillatory  convergence.  Note  that  solid  and  dotted  curves  are 
shown  for  each  order.  The  solid  curves  are  obtained  from  Equation  1,  as  shown,  with  the 

finite-difference  approximation  accurate  to  0(Ax2).  The  dotted  curves  are  obtained  from 

Equation  1  with  a  finite-difference  approximation  accurate  to  D(Ax^),  as  shown  in  the 
Appendix.  The  two  methods  show  nearly  identical  results  at  larger  order.  This  eliminates 
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the  possibility  of  the  finite-difference  approximation  being  the  cause  of  the  slow 
convergence.  At  this  point,  we  do  not  know  the  cause  of  the  slow  convergence.  In  the 
h/d  =10  column.  Table  1  again  compares  the  result  with  previous  results  from  References 
1  and  2.  The  agreement  is  excellent  for  TE  polarization.  For  TM  polarization,  the  previous 
results  are  not  reliable  as  discussed  in  the  last  paragraph.  There  is  no  other  exact  result  that 
can  be  used  to  check  for  this  deep  a  grating. 


Number  of  Orders 


FIGURE  4.  Relative  TM-Polarized  Diffraction  Intensity  of  the 

Reflecting  0,-1,  and  -2  Orders  Versus  Number  of  Orders.  The 
grating  profile  is  sinusoidal,  and  the  grating  material  is  an 

absorbing  metal  with  permittivity  (-48.91,  4.2).  The  grating 
height  is  h  =  llX,  period  d  =  1.7A,  angle  of  incidence  30  degrees, 
and  the  number  of  profile  sublayers  n^=  50.  This  case  is  identical 
to  Figure  2  except  the  grating  height  is  10  times  larger.  The  solid 
and  dotted  lines,  which  are  superimposed  for  the  larger  number  of 

orders,  are  for  finite-difference  approximations  accurate  to  O(Ax^) 
and  0(Ax‘^),  respectively. 
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Number  of  Orders 

FIGURE  5.  Relative  TE-Polarized  Diffraction  Intensity  of  the 

Reflecting  0,-1,  and  -2  Orders  Versus  Number  of  Orders.  This 
case  is  identical  to  Figure  4  except  the  incident  beam  is  TE  polarized. 


For  the  deep  grating  considered  in  Figures  4  and  5,  the  convergence  is  much  better  if 
we  change  the  permittivity  from  =  (-48.91,4.2)  to  =  (2.25,0)  to  represent  a 
nonabsorbing  dielectric.  This  grating  is  also  discussed  in  References  1  and  2.  For  this 
case.  Figures  6  and  7  show  that  the  convergence  is  very  good,  in  marked  contrast  to  the 
metallic  case.  For  this  dielectric  grating.  Table  2  compares  the  results  with  those  of 
References  1  and  2,  and  the  agreement  is  excellent  all  around. 
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Number  of  Orders 

FIGURE  6.  Relative  TM-Polarized  Diffraction  Intensity  for  the 

Transmitting  0,-1,  and  -2  Orders  Versus  Number  of  Orders.  The 
grating  profile  is  sinusoidal  and  the  material  is  a  nonabsorbing 
dielectric  with  permittivity  (2.25, 0).  The  grating  height  is  h  =  \l?i, 
period  d  =  \.lX,  angle  of  incidence  is  30  degrees,  and  the  number 
of  profile  sublayers  =  50  This  case  is  identical  to  Figure  4  except 
the  substrate  material  is  a  dielectric,  and  the  orders  are  transmitting. 


Number  of  Orders 

FIGURE  7.  Relative  TE-Polarized  Diffraction  Intensity  for  the 

Transmitting  0,-1,  and  -2  Orders  Versus  Number  of  Orders.  This 
case  is  identical  to  Figure  6  except  the  incident  beam  is  TE  polarized. 
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TABLE  2.  Transmission  Diffraction  Efficiencies 
Nonabsorbing  Dielectric  Sinusoidal  Grating. 


yd=  10 

Order 

5  point 
r-space 
average 

Our 

previous 

fc-space 

Li’s 

TE 

polarization 

-2 

0.4940 

0.4924 

0.4974 

-1 

0.4109 

0.4113 

0.4111 

0 

I>Td:g6WM 

TM 

polarization 

-2 

0.1867 

0.1905 

0.1869 

-1 

0.2744 

0.2812 

0.2757 

0 

0.5283 

0.5185 

0.5261 

In  conclusion,  we  described  a  real-space  modal  expansion  R-matrix  method,  and  we 
investigated  the  convergence  versus  number  of  diffracted  orders.  The  method  was  still 
stable  and  works  well  for  both  TM  and  TE  polarization  with  fast  convergence  for  absorbing 

metal  gratings  with  h  =  1.7X.  (h/d  -  1).  For  the  same  absorbing  metallic  grating  10  times 

deeper  Qi  =  \1X),  the  method  still  stable  for  both  TM  and  TE  polarization,  but  the 
convergence  appears  to  be  quite  slow  with  a  damped  oscillatory  behavior.  However,  when 
the  grating  material  is  replaced  by  a  dielectric,  the  convergence  becomes  quite  rapid.  We 
hope  that  future  studies  can  identify  the  cause  of  slow  convergence  for  the  deep  metallic 
grating. 
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Appendix 

In  deriving  Equation  1,  finite-difference  approximations  of  the  jc-derivatives  were 
made.  Consider  Ae  right-hand  side  of  Equation  la  where  the  centered-derivative 
approximation 


d(\dHy^ 

dx  j 


1  \Hy{x  +  ^x,z)-Hy{x,z)  Hy{x-  ^x,z)-  Hy{x,z) 

(Ax)^[  e{x  + i!al2,Zf.)  e(x- Ax/2, J 


(A-1) 


has  been  used  which  is  accurate  to  O(Ax^).  This  expression  is  obtained  by  writing  the 
formulas 


F(x  +  Ax/2)«F(x)-HyF'(x)  +  |(^yl  F"(x)  +  |fy  |  F'"(x) 


(A-2) 


Ay  1  f  Ay  1  f 

F{x-Ax/2)^F(x)-—F'(x)  +  -\^—j 


(A-3) 


and  solving  for  F'(x)  yields 


F'(x)  =  ^(^  +  ^/2)-F(x-Ax/2)  _  F"'ix)1^2 


Ax 


3  8 


(A-4) 


This  expression  is  clearly  accurate  to  order  O(Ax^).  In  Equation  A-4  we  may  replace  Ax 
by  3Ax  which  yields 


F'(x)« 


F(x-h3Ax/2)-F(x-3Ax/2) 

3Ax 


F^^^(x)  9 
3  8 


Ax^ 


(A-5) 


We  may  now  combine  Equations  A-4  and  A-5  such  that  the  O(Ax^)  term  vanishes  and  this 

yields  an  expression  for  F'(x),  which  is  accurate  to  order  OiAx'^).  This  0(Ax'')  formula 
may  be  applied  twice  to  evaluate  the  left-hand  side  of  Equation  A-1.  First  we  let 
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F(x)^Hyix)  and  this  yields  an  expression  F{(x)~dHy/dx.  Applying  the  formula  a 
second  time  with  F(x)  =  (l/e)F{{x)  yields  a  fourth-order  accurate  finite-difference  result 
F^ix)  =  [{H  e)F({x)\'  ~  {d!  dx)[(\l  £){dHy  /  dx)]  for  the  left-hand  side  of  Equation  A-1.  An 

analogous  procedure  can  be  applied  to  the  right-hand  side  of  Equation  Ic.  We  have  used 
both  die  second  and  fourth  order  formulas  in  our  numerical  calculations.  Unfortunately, 
however,  no  significant  improvement  in  convergence  was  seen  for  the  deep  grating 
calculations  considered  here. 
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